library(tidyverse)
library(gtsummary)
library(flextable)
library(corrplot)
library(pacman)
library(devtools)
library(qgcomp)
library(forcats)
library(cowplot)
library(broom)
library(bkmr)
library(qgcompint)
library(kohonen)
library(sommix)
library(prettydoc)
pacman::p_load(future, future.apply, reshape2, MASS, gWQS)
Distribution of sociodemographic characteristics and Child Behavior Checklist (CBCL) outcomes among Black mother-child pairs in metropolitan Atlanta, Georgia.
# Select variables
table_1 <- mix %>% dplyr::select(`Maternal Age`, `Maternal Education`, `Marital Status`, Parity,`Pre-pregnancy BMI`,`Substance Use`, `Insurance Type`, `Income to Poverty Ratio`, `Infant Sex`, `Mean Child Age`, Internalizing, Externalizing)
# Set theme
theme_gtsummary_journal(journal = "lancet")
# Table
table_1 <- table_1 %>%
tbl_summary(type = all_dichotomous() ~ "categorical", statistic = all_continuous() ~ "{mean} ({sd})", digits = all_continuous() ~ 2) %>%
modify_header(label = "**Demographic**") %>%
bold_labels()
table_1
| Demographic | N = 1591 |
|---|---|
| Maternal Age | 24·97 (4·46) |
| Maternal Education | |
| Â Â Â Â Less than some college | 83 (52%) |
| Â Â Â Â College degree or higher | 76 (48%) |
| Marital Status | |
| Â Â Â Â Never married | 75 (47%) |
| Â Â Â Â Married | 84 (53%) |
| Parity | |
| Â Â Â Â No prior births | 65 (41%) |
| Â Â Â Â 1 or more prior births | 94 (59%) |
| Pre-pregnancy BMI | |
| Â Â Â Â Underweight or normal | 64 (40%) |
| Â Â Â Â Overweight | 29 (18%) |
| Â Â Â Â Obese | 66 (42%) |
| Substance Use | |
| Â Â Â Â No | 90 (57%) |
| Â Â Â Â Yes | 69 (43%) |
| Insurance Type | |
| Â Â Â Â Low-income Medicaid | 69 (43%) |
| Â Â Â Â Pregnancy Medicaid | 58 (36%) |
| Â Â Â Â Private | 32 (20%) |
| Income to Poverty Ratio | |
| Â Â Â Â <100% | 76 (48%) |
| Â Â Â Â 100% - 199% | 52 (33%) |
| Â Â Â Â >200% | 31 (19%) |
| Infant Sex | |
| Â Â Â Â Male | 81 (51%) |
| Â Â Â Â Female | 78 (49%) |
| Mean Child Age | 2·82 (0·64) |
| Internalizing | 7·10 (5·29) |
| Externalizing | 11·37 (7·62) |
| 1 Mean (SD); n (%) | |
Distributions of per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in serum samples obtained during the first trimester, among Black pregnant women in metropolitan Atlanta, Georgia.
# Select variables
table_s1 <- mix %>% dplyr::select(PFHXS, PFOS, PFOA, PFNA, PBDE47, PBDE99, PBDE100)
table_s1_PFAS <- atl_chemicals %>% dplyr::select(PFHXS, PFOS, PFOA, PFNA) %>% drop_na(PFHXS, PFOS, PFOA, PFNA)
table_s1_PBDE <- atl_chemicals %>% dplyr::select(PBDE47, PBDE99, PBDE100) %>% drop_na(PBDE47, PBDE99, PBDE100)
# Build function
dist.func <- function(x) {
c(
geomean = round(exp(mean(log(x), na.rm = TRUE)), 2),
geosd = round(exp(sd(log(x), na.rm = TRUE)), 2),
min = round(min(x, na.rm = TRUE), 2),
max = round(max(x, na.rm = TRUE), 2),
median = round(median(x, na.rm = TRUE), 2)
)
}
# Apply
table_s1 <- as.data.frame(do.call(rbind, lapply(table_s1, dist.func)))
table_s1 <- table_s1 %>%
mutate(sample = "Analytic")
table_s1_PFAS <- as.data.frame(do.call(rbind, lapply(table_s1_PFAS, dist.func)))
table_s1_PFAS <- table_s1_PFAS %>%
mutate(sample = "Full")
table_s1_PBDE <- as.data.frame(do.call(rbind, lapply(table_s1_PBDE, dist.func)))
table_s1_PBDE <- table_s1_PBDE %>%
mutate(sample = "Full")
table_s1 <- rbind(table_s1_PFAS, table_s1, table_s1_PBDE)
# Table
table_s1 <- table_s1 %>%
mutate(Chemical = rownames(table_s1)) %>%
mutate("Geometric Mean (SD)" = paste0(geomean, " (", geosd, ")")) %>%
dplyr::select(Chemical, sample, `Geometric Mean (SD)`, min, median, max) %>%
as_grouped_data("sample") %>%
flextable() %>%
set_header_labels(
geomean = "Geometric Mean",
geosd = "Geometric SD",
min = "Minimum",
max = "Maximum",
median = "Median") %>%
font(fontname = "Tahoma") %>%
autofit() %>%
theme_alafoli()
table_s1
sample | Chemical | Geometric Mean (SD) | Minimum | Median | Maximum |
|---|---|---|---|---|---|
Full | |||||
PFHXS | 1.16 (2.06) | 0.08 | 1.24 | 6.17 | |
PFOS | 1.89 (2.5) | 0.01 | 2.12 | 12.42 | |
PFOA | 0.62 (2.36) | 0.01 | 0.70 | 4.42 | |
PFNA | 0.26 (2.35) | 0.01 | 0.30 | 2.27 | |
Analytic | |||||
PFHXS1 | 1.07 (2.07) | 0.14 | 1.27 | 4.50 | |
PFOS1 | 2.37 (1.97) | 0.14 | 2.53 | 12.42 | |
PFOA1 | 0.8 (2.12) | 0.03 | 0.89 | 4.42 | |
PFNA1 | 0.33 (2.04) | 0.01 | 0.34 | 2.27 | |
PBDE47 | 92.7 (1.98) | 20.65 | 87.48 | 441.86 | |
PBDE99 | 23.07 (2.21) | 7.81 | 23.26 | 308.29 | |
PBDE100 | 14.6 (2.85) | 2.81 | 17.28 | 111.47 | |
Full | |||||
PBDE471 | 89.84 (2.08) | 15.12 | 83.86 | 1,382.08 | |
PBDE991 | 23.1 (2.35) | 7.81 | 22.11 | 341.98 | |
PBDE1001 | 13.32 (2.99) | 2.81 | 15.44 | 289.20 |
Pearson correlation for the association between natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in first trimester serum samples, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.
# Correlation matrix and plot
mix_cor <- mix %>% dplyr::select(PFHXS, PFOS, PFOA, PFNA, PBDE47, PBDE99, PBDE100, Internalizing, Externalizing)
cor <- cor(mix_cor, method = "pearson")
figure_s3 <- corrplot(cor, method = "color", tl.col = "grey25", tl.cex = 6, cl.cex = 6, mar=c(0,0,1,0)) %>%
corrRect(c(1,5,8, 9))
# Total Sample
# Assign predictors, outcomes
predictors <- c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")
outcomes <- c("Internalizing", "Externalizing")
# Loop regression
models <- list()
for (i in 1:length(outcomes)) {
for (j in 1:length(predictors)) {
model_name <- paste0(outcomes[i], " ~ ", predictors[j])
formula <- as.formula(paste(outcomes[i], " ~ ", predictors[j], " + `Mean Child Age` + `Infant Sex` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`"))
model <- lm(formula = formula, data = mix)
models[[model_name]] <- tidy(model, conf.int = TRUE) %>%
mutate(outcome = outcomes[i])
}
}
chemicals_adj <- do.call(rbind, models)
# Format
reg <- chemicals_adj %>%
filter(term %in% c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")) %>%
dplyr::select(outcome, term, estimate, conf.low, conf.high) %>%
mutate(sex = "All")
### By sex
# Male
## Filter for males
mix_m <- mix %>% filter(`Infant Sex` == "Male")
models <- list()
for (i in 1:length(outcomes)) {
for (j in 1:length(predictors)) {
model_name <- paste0(outcomes[i], " ~ ", predictors[j])
formula <- as.formula(paste(outcomes[i], " ~ ", predictors[j], " + `Mean Child Age` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`"))
model <- lm(formula = formula, data = mix_m)
models[[model_name]] <- tidy(model, conf.int = TRUE) %>%
mutate(outcome = outcomes[i])
}
}
chemicals_adj <- do.call(rbind, models)
reg_male <- chemicals_adj %>%
filter(term %in% c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")) %>%
dplyr::select(outcome, term, estimate, conf.low, conf.high) %>%
mutate(sex = "Male")
# Female
## Filter for females
mix_f <- mix %>% filter(`Infant Sex` == "Female")
models <- list()
for (i in 1:length(outcomes)) {
for (j in 1:length(predictors)) {
model_name <- paste0(outcomes[i], " ~ ", predictors[j])
formula <- as.formula(paste(outcomes[i], " ~ ", predictors[j], " + `Mean Child Age` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`"))
model <- lm(formula = formula, data = mix_f)
models[[model_name]] <- tidy(model, conf.int = TRUE) %>%
mutate(outcome = outcomes[i])
}
}
chemicals_adj <- do.call(rbind, models)
reg_female <- chemicals_adj %>%
filter(term %in% c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")) %>%
dplyr::select(outcome, term, estimate, conf.low, conf.high) %>%
mutate(sex = "Female")
Associations between natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia. PFAS in ng.mL and PBDEs in pg/mL.
Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use
# Format
reg$outcome <- factor(reg$outcome, levels = c("Internalizing", "Externalizing"))
reg$term <- factor(reg$term, levels = c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))
chemical_label <- c("PFHXS", "PFOS", "PFOA", "PFNA", "BDE-47", "BDE-99", "BDE-100")
names(chemical_label) <- c("PFHXS.log", "PFOS.log","PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")
# Figure
figure_1 <- ggplot(reg, aes(x = estimate, y = term, color = outcome)) +
geom_pointrange(aes(xmin = conf.low, xmax = conf.high), size = 4, linewidth = 4) +
geom_vline(xintercept = 0, color = "grey25", linetype = "dashed", alpha = 0.5) +
facet_grid(~outcome, scales = "free") +
xlab("Beta (95% Confidence Interval)") +
ylab("Chemical") +
scale_x_continuous(n.breaks=4) +
scale_y_discrete(labels = chemical_label) +
scale_color_manual(values = c("grey25", "grey55")) +
coord_flip() +
theme_classic(base_family = "AppleGothic") +
theme(title = element_text(size = 50, color = "grey25"),
axis.title.x = element_text(size = 50, color = "grey25"),
axis.title.y = element_text(size = 50, color = "grey25"),
axis.text.x = element_text(size = 50, color = "grey25", angle = 45, vjust = 0.6),
axis.text.y = element_text(size = 50, color = "grey25"),
strip.background = element_rect(fill="white"),
strip.text = element_text(colour = 'grey25', face = "bold", size = 50),
legend.position = "none",
strip.background.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
figure_1
Linear regression associations (overall, child sex specific) between natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.
Total sample models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use. Sex specific models are adjusted for mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.
# Combine regression results
reg <- rbind(reg, reg_female, reg_male)
# Table
table_s2 <- reg %>%
mutate_if(is.numeric, round, digits = 2) %>%
mutate(Beta_CI = paste0(estimate, " (", conf.low, ", ", conf.high, ")")) %>%
mutate(term = recode(term,
"PFHXS.log" = "PFHXS",
"PFOA.log" = "PFOA",
"PFNA.log" = "PFNA",
"PFOS.log" = "PFOS",
"PBDE47.log" = "BDE-47",
"PBDE99.log" = "BDE-99",
"PBDE100.log" = "BDE-100")) %>%
dplyr::select(term, outcome, sex, Beta_CI) %>%
pivot_wider(names_from = outcome, values_from = (Beta_CI)) %>%
group_by(term, sex) %>%
summarise(
Internalizing = Internalizing,
Externalizing = Externalizing) %>%
as_grouped_data(groups = "term") %>%
flextable() %>%
add_body_row(values = c("","" ,"Beta (95% CI)", "Beta (95% CI)")) %>%
font(fontname = "Tahoma") %>%
set_header_labels(term = "Chemical", sex = "Sex") %>%
autofit() %>%
theme_alafoli()
table_s2
Chemical | Sex | Internalizing | Externalizing |
|---|---|---|---|
Beta (95% CI) | Beta (95% CI) | ||
PFHXS | |||
All | -0.17 (-0.33, 0) | -0.21 (-0.36, -0.05) | |
Female | -0.12 (-0.35, 0.11) | -0.17 (-0.37, 0.02) | |
Male | -0.17 (-0.42, 0.07) | -0.22 (-0.48, 0.03) | |
PFOS | |||
All | -0.1 (-0.26, 0.06) | -0.21 (-0.36, -0.05) | |
Female | -0.01 (-0.25, 0.23) | -0.13 (-0.33, 0.06) | |
Male | -0.16 (-0.39, 0.06) | -0.24 (-0.48, -0.01) | |
PFOA | |||
All | -0.12 (-0.28, 0.04) | -0.17 (-0.33, -0.01) | |
Female | -0.02 (-0.36, 0.32) | -0.13 (-0.42, 0.15) | |
Male | -0.18 (-0.39, 0.02) | -0.21 (-0.43, 0) | |
PFNA | |||
All | -0.15 (-0.32, 0.01) | -0.21 (-0.36, -0.05) | |
Female | -0.04 (-0.29, 0.2) | -0.06 (-0.27, 0.15) | |
Male | -0.23 (-0.46, -0.01) | -0.31 (-0.55, -0.07) | |
BDE-47 | |||
All | 0.11 (-0.06, 0.28) | 0.16 (-0.01, 0.32) | |
Female | 0.02 (-0.2, 0.25) | 0 (-0.19, 0.18) | |
Male | 0.22 (-0.04, 0.48) | 0.31 (0.04, 0.58) | |
BDE-99 | |||
All | 0.09 (-0.07, 0.25) | 0.15 (-0.01, 0.3) | |
Female | 0.1 (-0.1, 0.31) | 0.13 (-0.04, 0.3) | |
Male | 0.09 (-0.16, 0.34) | 0.14 (-0.12, 0.41) | |
BDE-100 | |||
All | 0.07 (-0.1, 0.23) | 0.15 (-0.01, 0.31) | |
Female | 0.03 (-0.2, 0.26) | 0.05 (-0.14, 0.25) | |
Male | 0.08 (-0.17, 0.34) | 0.17 (-0.1, 0.44) |
Linear regression associations between natural log transformed polybrominated diphenyl ethers (PBDEs), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia. Sensitivity analysis additionally adjusted for total lipids.
# Assign predictors, outcomes
predictors <- c("PBDE47.log", "PBDE99.log", "PBDE100.log")
outcomes <- c("Internalizing", "Externalizing")
# Loop regression
models <- list()
for (i in 1:length(outcomes)) {
for (j in 1:length(predictors)) {
model_name <- paste0(outcomes[i], " ~ ", predictors[j])
formula <- as.formula(paste(outcomes[i], " ~ ", predictors[j], " + `Mean Child Age` + `Infant Sex` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use` + total_lipids"))
model <- lm(formula = formula, data = mix)
models[[model_name]] <- tidy(model, conf.int = TRUE) %>%
mutate(N = nobs(model)) %>%
mutate(outcome = outcomes[i])
}
}
chemicals_adj_lipid <- do.call(rbind, models)
# Format
reg_lipid <- chemicals_adj_lipid %>%
filter(term %in% c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")) %>%
dplyr::select(outcome, N, term, estimate, conf.low, conf.high)
# Combine regression results
# Table
table_s3 <- reg_lipid %>%
mutate_if(is.numeric, round, digits = 2) %>%
mutate(`Beta CI` = paste0(estimate, " (", conf.low, ", ", conf.high, ")")) %>%
mutate(term = recode(term,
"PFHXS.log" = "PFHXS",
"PFOA.log" = "PFOA",
"PFNA.log" = "PFNA",
"PFOS.log" = "PFOS",
"PBDE47.log" = "BDE-47",
"PBDE99.log" = "BDE-99",
"PBDE100.log" = "BDE-100")) %>%
dplyr::select(term, N, outcome, `Beta CI`) %>%
pivot_wider(names_from = outcome, values_from = c("N", "Beta CI"),names_glue = "{outcome}_{.value}") %>%
dplyr::select(term, Internalizing_N, `Internalizing_Beta CI`, Externalizing_N, `Externalizing_Beta CI`) %>%
flextable() %>%
separate_header() %>%
font(fontname = "Tahoma") %>%
set_header_labels(term = "Chemical") %>%
autofit() %>%
theme_alafoli()
table_s3
term | Internalizing | Externalizing | ||
|---|---|---|---|---|
N | Beta CI | N | Beta CI | |
BDE-47 | 154 | 0.12 (-0.05, 0.29) | 154 | 0.16 (-0.01, 0.33) |
BDE-99 | 154 | 0.09 (-0.07, 0.25) | 154 | 0.14 (-0.02, 0.3) |
BDE-100 | 154 | 0.09 (-0.08, 0.25) | 154 | 0.17 (0, 0.33) |
Linear regression associations between natural log transformed per- and polyfluoroalkyl substances (PFAS) and polybrominated diphenyl ethers (PBDEs), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia. Sensitivity analysis removing preterm births.
# Removing preterm brith
mix_ptb <- mix %>% filter(ptb == 0)
# Assign predictors, outcomes
predictors <- c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")
outcomes <- c("Internalizing", "Externalizing")
# Loop regression
models <- list()
for (i in 1:length(outcomes)) {
for (j in 1:length(predictors)) {
model_name <- paste0(outcomes[i], " ~ ", predictors[j])
formula <- as.formula(paste(outcomes[i], " ~ ", predictors[j], " + `Mean Child Age` + `Infant Sex` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`"))
model <- lm(formula = formula, data = mix_ptb)
models[[model_name]] <- tidy(model, conf.int = TRUE) %>%
mutate(outcome = outcomes[i])
}
}
chemicals_adj_ptb <- do.call(rbind, models)
# Format
reg_ptb <- chemicals_adj_ptb %>%
filter(term %in% c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")) %>%
dplyr::select(outcome, term, estimate, conf.low, conf.high)
# Combine regression results
# Table
table_s4 <- reg_ptb %>%
mutate_if(is.numeric, round, digits = 2) %>%
mutate(Beta_CI = paste0(estimate, " (", conf.low, ", ", conf.high, ")")) %>%
mutate(term = recode(term,
"PFHXS.log" = "PFHXS",
"PFOA.log" = "PFOA",
"PFNA.log" = "PFNA",
"PFOS.log" = "PFOS",
"PBDE47.log" = "BDE-47",
"PBDE99.log" = "BDE-99",
"PBDE100.log" = "BDE-100")) %>%
dplyr::select(term, outcome, Beta_CI) %>%
pivot_wider(names_from = outcome, values_from = (Beta_CI)) %>%
group_by(term) %>%
summarise(
Internalizing = Internalizing,
Externalizing = Externalizing) %>%
flextable() %>%
add_body_row(values = c("","Beta (95% CI)", "Beta (95% CI)")) %>%
font(fontname = "Tahoma") %>%
set_header_labels(term = "Chemical") %>%
autofit() %>%
theme_alafoli()
table_s4
Chemical | Internalizing | Externalizing |
|---|---|---|
Beta (95% CI) | Beta (95% CI) | |
BDE-100 | 0.05 (-0.13, 0.22) | 0.13 (-0.04, 0.31) |
BDE-47 | 0.17 (-0.02, 0.35) | 0.16 (-0.03, 0.35) |
BDE-99 | 0.13 (-0.04, 0.3) | 0.15 (-0.03, 0.32) |
PFHXS | -0.18 (-0.35, -0.01) | -0.24 (-0.41, -0.07) |
PFNA | -0.16 (-0.33, 0.01) | -0.22 (-0.4, -0.05) |
PFOA | -0.11 (-0.29, 0.06) | -0.2 (-0.37, -0.03) |
PFOS | -0.12 (-0.29, 0.05) | -0.24 (-0.41, -0.07) |
# Define exposure, outcome, dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, `Pre-pregnancy BMI`, `Substance Use`)
exposure.names<-c("PFHXS.log", "PFOS.log", "PFOA.log","PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")
outcomes <- c("Internalizing", "Externalizing")
# Build function for qgcomp
qgcompdata <- function(outcome) {
model <- qgcomp.glm.noboot(as.formula(paste(outcome, "~", "PFHXS.log + PFOS.log + PFOA.log + PFNA.log + PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` + +`Substance Use`")),
expnms = exposure.names,
dat = mixture,
family = gaussian(),
q = 4)
coefs <- data.frame(
estimate = round(model$psi,2),
lower_ci = round(model$ci.coef[-1, 1],2),
upper_ci = round(model$ci.coef[-1, 2],2),
outcome = outcome,
row.names = NULL
)
adj_df <- data.frame()
adj_df <- rbind(adj_df, coefs)
print(adj_df)
}
# Run
results_all_int <- qgcompdata(outcome = "Internalizing")
results_all_ext <- qgcompdata(outcome = "Externalizing")
# Combine, format
qgcomp_all_df <- rbind(results_all_int, results_all_ext)
qgcomp_all_df <- qgcomp_all_df %>% mutate(group = "Total Mixture")
##### PFAS
# Define exposure, outcome dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, `Pre-pregnancy BMI`, `Insurance Type`, `Substance Use`, `Marital Status`)
outcomes <- c("Internalizing", "Externalizing")
exposure.names<-c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log")
# Function for qgcomp
qgcompdata <- function(outcome) {
model <- qgcomp.glm.noboot(as.formula(paste(outcome, "~", "PFHXS.log + PFOS.log + PFOA.log + PFNA.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` + `Insurance Type` +`Substance Use` +`Marital Status`")),
expnms = exposure.names,
dat = mixture,
family = gaussian(),
q = 4)
coefs <- data.frame(
estimate = round(model$psi,2),
lower_ci = round(model$ci.coef[-1, 1],2),
upper_ci = round(model$ci.coef[-1, 2],2),
outcome = outcome,
row.names = NULL
)
adj_df <- data.frame()
PFAS_df <- rbind(adj_df, coefs)
print(PFAS_df)
}
# Run
results_pfas_int <- qgcompdata(outcome = "Internalizing")
results_pfas_ext <- qgcompdata(outcome = "Externalizing")
qgcomp_pfas_df <- rbind(results_pfas_int, results_pfas_ext)
#### PBDE
# Define exposure, outcome, dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, `Pre-pregnancy BMI`, `Insurance Type`, `Substance Use`, `Marital Status`)
exposure.names<-c("PBDE47.log", "PBDE99.log", "PBDE100.log")
outcomes <- c("Internalizing", "Externalizing")
qgcompdata <- function(outcome) {
model <- qgcomp.glm.noboot(as.formula(paste(outcome, "~", "PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` + `Insurance Type` + +`Substance Use` +`Marital Status`")),
expnms = exposure.names,
dat = mixture,
family = gaussian(),
q = 4)
coefs <- data.frame(
estimate = round(model$psi,2),
lower_ci = round(model$ci.coef[-1, 1],2),
upper_ci = round(model$ci.coef[-1, 2],2),
outcome = outcome,
row.names = NULL
)
adj_df <- data.frame()
adj_df <- rbind(adj_df, coefs)
print(adj_df)
}
# Run
results_pbde_int <- qgcompdata(outcome = "Internalizing")
results_pbde_ext <- qgcompdata(outcome = "Externalizing")
# Combine results
qgcomp_pbde_df <- rbind(results_pbde_int, results_pbde_ext)
qgcomp_chemicals_df <- rbind(qgcomp_pbde_df, qgcomp_pfas_df)
qgcomp_chemicals_df <- qgcomp_chemicals_df %>% mutate(group = c("PBDE", "PBDE", "PFAS", "PFAS"))
Quantile g-computation cumulative effect (95% confidence intervals) and weights reflecting the associations between a mixture of natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.
Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use. multiply by negative 1
##### TOTAL SAMPLE
# Forest plot
qgcomp_all_df$outcome <- factor(qgcomp_all_df$outcome, levels = c("Internalizing", "Externalizing"))
forest_all_qgcomp <- ggplot(qgcomp_all_df, aes(estimate, outcome)) +
geom_pointrange(aes(xmin = lower_ci, xmax = upper_ci),
size = 2, linewidth = 2) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey25") +
coord_flip() +
scale_x_continuous(limits = c(-0.6, 0.4), breaks = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4)) +
ggtitle("Total Sample") +
theme_classic(base_family = "AppleGothic") +
labs(x = "Beta (95 % Confidence Interval)", y = "CBCL Outcome") +
theme(text = element_text(size = 25, color = "grey25", family = "Helvetica"),
legend.title = element_text(color = "grey25"),
axis.text = element_text(color = "grey25", size = 25),
plot.title = element_text(hjust = 0.5))
# Weights
# Internalizing
exposure.names<-c("PFHXS.log", "PFOS.log", "PFOA.log","PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")
model <- qgcomp.glm.noboot(Internalizing ~ PFHXS.log + PFOS.log + PFOA.log + PFNA.log + PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
expnms = exposure.names,
dat = mixture,
family = gaussian(),
q = 4)
# Extract weights
pos <- as.data.frame(model$pos.weights)
# Format
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
mutate(chemical = c("BDE-47", "PFOS", "BDE-99", "PFOA"), outcome = "Internalizing", Direction = "Positive")
neg <- as.data.frame(model$neg.weights)
neg <- neg %>% rename(weight = `model$neg.weights`) %>%
mutate(chemical = c("PFHXS", "PFNA", "BDE-100"), outcome = "Internalizing", Direction = "Negative")
neg$weight <- neg$weight*-1
int_all_weights_qgcomp <- rbind(pos, neg)
# Externalizing
model <- qgcomp.glm.noboot(Externalizing ~ PFHXS.log + PFOS.log + PFOA.log + PFNA.log + PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
expnms = exposure.names,
dat = mixture,
family = gaussian(),
q = 4)
# Extract weights
pos <- as.data.frame(model$pos.weights)
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
mutate(chemical = c("BDE-47", "PFOA", "BDE-100", "BDE-99", "PFOS"), outcome = "Externalizing", Direction = "Positive")
neg <- as.data.frame(model$neg.weights)
neg <- neg %>% rename(weight = `model$neg.weights`) %>%
mutate(chemical = c("PFHXS", "PFNA"), outcome = "Externalizing", Direction = "Negative")
neg$weight <- neg$weight*-1
ext_all_weights_qgcomp <- rbind(pos, neg)
# Plot
int_all_weights_qgcomp <- ggplot(int_all_weights_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
geom_col() +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Internalizing") +
scale_fill_manual(values = c("salmon3", "steelblue3")) +
scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
theme_classic(base_family = "AppleGothic") +
theme(text = element_text(size = 25, color = "grey25"),
legend.title = element_text(color = "grey25"),
axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
axis.text.y = element_text(color = "grey25", size = 25),
legend.position = "none",
plot.title = element_text(hjust = 0.5))
ext_all_weights_qgcomp <- ggplot(ext_all_weights_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
geom_col() +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Externalizing") +
scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
scale_fill_manual(values = c("salmon3", "steelblue3")) +
theme_classic(base_family = "AppleGothic") +
theme(text = element_text(size = 25, color = "grey25"),
legend.title = element_text(color = "grey25"),
axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
axis.text.y = element_text(color = "grey25", size = 25),
legend.position = "bottom",
plot.title = element_text(hjust = 0.5),
legend.key.size = unit(3, "cm"),
legend.text = element_text(size = 25))
##### CHEMCIAL SPECIFIC
# Forest plot
qgcomp_pfas_df$outcome <- factor(qgcomp_pfas_df$outcome, levels = c("Internalizing", "Externalizing"))
forest_pfas_qgcomp <- ggplot(qgcomp_pfas_df, aes(estimate, outcome)) +
geom_pointrange(aes(xmin = lower_ci, xmax = upper_ci),
size = 2, linewidth = 2, color = "black", position = position_dodge(width = 0.25)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey25") +
scale_x_continuous(limits = c(-0.6, 0.4), breaks = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4)) +
ggtitle("PFAS") +
coord_flip() +
theme_classic(base_family = "AppleGothic") +
labs(x = "Beta (95 % Confidence Interval)", y = "CBCL Outcome", linetype = "Chemical Class") +
theme(text = element_text(size = 25, color = "grey25", family = "Helvetica"),
legend.title = element_text(color = "grey25"),
axis.text = element_text(color = "grey25", size = 25),
plot.title = element_text(hjust = 0.5))
forest_pbde_qgcomp <- ggplot(qgcomp_pbde_df, aes(estimate, outcome)) +
geom_pointrange(aes(xmin = lower_ci, xmax = upper_ci),
size = 2, linewidth = 2, color = "black", position = position_dodge(width = 0.25)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey25") +
scale_x_continuous(limits = c(-0.6, 0.4), breaks = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4)) +
ggtitle("PBDE") +
coord_flip() +
theme_classic(base_family = "AppleGothic") +
labs(x = "Beta (95 % Confidence Interval)", y = "CBCL Outcome", linetype = "Chemical Class") +
theme(text = element_text(size = 25, color = "grey25", family = "Helvetica"),
legend.title = element_text(color = "grey25"),
axis.text = element_text(color = "grey25", size = 25),
plot.title = element_text(hjust = 0.5))
# Weights
# PFAS
exposure.names<-c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log")
model <- qgcomp.glm.noboot(Internalizing ~ PFHXS.log + PFOS.log + PFOA.log + PFNA.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
expnms = exposure.names,
dat = mixture,
family = gaussian(),
q = 4)
# Extract weights
pos <- as.data.frame(model$pos.weights)
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
mutate(chemical = c("PFOS", "PFOA"), outcome = "Internalizing", Direction = "Positive")
neg <- as.data.frame(model$neg.weights)
neg <- neg %>% rename(weight = `model$neg.weights`) %>%
mutate(chemical = c("PFNA", "PFHXS"), outcome = "Internalizing", Direction = "Negative")
neg$weight <- neg$weight*-1
weights_int_pfas_qgcomp <- rbind(pos, neg)
# Externalizing
model <- qgcomp.glm.noboot(Externalizing ~ PFHXS.log + PFOS.log + PFOA.log + PFNA.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
expnms = exposure.names,
dat = mixture,
family = gaussian(),
q = 4)
# Extract weights
pos <- as.data.frame(model$pos.weights)
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
mutate(chemical = c("PFOA", "PFOS"), outcome = "Externalizing", Direction = "Positive")
neg <- as.data.frame(model$neg.weights)
neg <- neg %>% rename(weight = `model$neg.weights`) %>%
mutate(chemical = c("PFNA", "PFHXS"), outcome = "Externalizing", Direction = "Negative")
neg$weight <- neg$weight*-1
weights_ext_pfas_qgcomp <- rbind(pos, neg)
# PBDE
exposure.names<-c("PBDE47.log", "PBDE99.log", "PBDE100.log")
model <- qgcomp.glm.noboot(Internalizing ~ PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
expnms = exposure.names,
dat = mixture,
family = gaussian(),
q = 4)
pos <- as.data.frame(model$pos.weights)
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
mutate(chemical = c("BDE-47", "BDE-99"), outcome = "Internalizing", Direction = "Positive")
neg <- as.data.frame(model$neg.weights)
neg <- neg %>% rename(weight = `model$neg.weights`) %>%
mutate(chemical = c("BDE-100"), outcome = "Internalizing", Direction = "Negative")
neg$weight <- neg$weight*-1
# Extract weights
weights_int_pbde_qgcomp <- rbind(pos, neg)
# Externalizing
model <- qgcomp.glm.noboot(Externalizing ~ PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
expnms = exposure.names,
dat = mixture,
family = gaussian(),
q = 4)
# Extract weights
pos <- as.data.frame(model$pos.weights)
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
mutate(chemical = c("BDE-47", "BDE-99", "BDE-100"), outcome = "Externalizing", Direction = "Positive")
weights_ext_pbde_qgcomp <- pos
# Plot
int_pfas_weights_qgcomp <- ggplot(weights_int_pfas_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
geom_col() +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Internalizing") +
scale_fill_manual(values = c("salmon3", "steelblue3")) +
scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
theme_classic(base_family = "AppleGothic") +
theme(text = element_text(size = 20, color = "grey25"),
legend.title = element_text(color = "grey25"),
axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
axis.text.y = element_text(color = "grey25", size = 25),
legend.position = "none",
plot.title = element_text(hjust = 0.5))
ext_pfas_weights_qgcomp <- ggplot(weights_ext_pfas_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
geom_col() +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Externalizing") +
scale_fill_manual(values = c("salmon3", "steelblue3")) +
scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
theme_classic(base_family = "AppleGothic") +
theme(text = element_text(size = 25, color = "grey25"),
legend.title = element_text(color = "grey25"),
axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
axis.text.y = element_text(color = "grey25", size = 25),
legend.position = "bottom",
plot.title = element_text(hjust = 0.5),
legend.key.size = unit(3, "cm"),
legend.text = element_text(size = 25))
int_pbde_weights_qgcomp <- ggplot(weights_int_pbde_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
geom_col() +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Internalizing") +
scale_fill_manual(values = c("salmon3", "steelblue3")) +
theme_classic(base_family = "AppleGothic") +
scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
theme(text = element_text(size = 25, color = "grey25"),
legend.title = element_text(color = "grey25"),
axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
axis.text.y = element_text(color = "grey25", size = 25),
legend.position = "none",
plot.title = element_text(hjust = 0.5))
ext_pbde_weights_qgcomp <- ggplot(weights_ext_pbde_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
geom_col() +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Externalizing") +
scale_fill_manual(values = c("steelblue3")) +
theme_classic(base_family = "AppleGothic") +
scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
theme(text = element_text(size = 25, color = "grey25"),
legend.title = element_text(color = "grey25"),
axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
axis.text.y = element_text(color = "grey25", size = 25),
legend.position = "bottom",
plot.title = element_text(hjust = 0.5),
legend.key.size = unit(3, "cm"),
legend.text = element_text(size = 25))
figure_3 <- plot_grid(forest_all_qgcomp, forest_pfas_qgcomp, forest_pbde_qgcomp, int_all_weights_qgcomp,
int_pfas_weights_qgcomp, int_pbde_weights_qgcomp,ext_all_weights_qgcomp, ext_pfas_weights_qgcomp,
ext_pbde_weights_qgcomp, ncol = 3, nrow = 3, align = "v")
figure_3
Association between the total mixture of per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in first trimester serum samples, a mixture of PFAS alone, and a mixture of PBDEs alone, estimated using quantile g-computation, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.
Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.
# Combine results
table_s6 <- rbind(qgcomp_all_df, qgcomp_chemicals_df)
# Table
table_s6 <- table_s6 %>%
mutate_if(is.numeric, round, digits = 2) %>%
mutate("Beta (CI)" = paste0(estimate, " (", lower_ci, ", ", upper_ci, ")")) %>%
dplyr::select(group, outcome, `Beta (CI)`) %>%
pivot_wider(names_from = outcome, values_from = `Beta (CI)`) %>%
flextable() %>%
font(fontname = "Tahoma") %>%
add_body_row(values = c("","Beta (95% CI)", "Beta (95% CI)")) %>%
merge_v(j = "group") %>%
set_header_labels(
group = "") %>%
autofit() %>%
theme_alafoli()
table_s6
Internalizing | Externalizing | |
|---|---|---|
Beta (95% CI) | Beta (95% CI) | |
Total Mixture | -0.04 (-0.28, 0.2) | -0.06 (-0.3, 0.18) |
PBDE | 0.13 (-0.03, 0.3) | 0.2 (0.04, 0.36) |
PFAS | -0.2 (-0.39, -0.01) | -0.28 (-0.46, -0.1) |
# Define dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi, `Substance Use`)
# Define mixture, covariates
df.exp <- mixture %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log)
df.cov <- mixture %>% dplyr::select(`Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi, `Substance Use`)
df.cov <- sapply(df.cov, as.numeric)
# BKMR
# Internalizing
set.seed(123)
fit_int <- kmbayes(y = mixture$Internalizing,
Z = df.exp,
X = df.cov,
iter = 10000,
verbose = FALSE,
varsel = TRUE,
family='gaussian')
# Extract pips
pips_int <- ExtractPIPs(fit_int)
pips_int <- pips_int %>% mutate(outcome = "Internalizing", group = "Total Sample")
# Externalizing
set.seed(123)
fit_ext <- kmbayes(y = mixture$Externalizing,
Z = df.exp,
X = df.cov,
iter = 10000,
verbose = FALSE,
varsel = TRUE,
family='gaussian')
# Extract pips
pips_ext <- ExtractPIPs(fit_ext)
pips_ext <- pips_ext %>% mutate(outcome = "Externalizing", group = "Total Sample")
## PFAS
# Define dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi, `Substance Use`)
# Define mixture, covariates
df.exp <- mixture %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log)
df.cov <- mixture %>% dplyr::select(`Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi,`Substance Use`)
df.cov <- sapply(df.cov, as.numeric)
# BKMR
# Internalizing
set.seed(123)
fit_int_pfas <- kmbayes(y = mixture$Internalizing,
Z = df.exp,
X = df.cov,
iter = 10000,
verbose = FALSE,
varsel = TRUE,
family='gaussian')
# Extract pips
pips_int_pfas <- ExtractPIPs(fit_int_pfas)
pips_int_pfas <- pips_int_pfas %>% mutate(outcome = "Internalizing", group = "PFAS")
# Externalizing
set.seed(123)
fit_ext_pfas <- kmbayes(y = mixture$Externalizing,
Z = df.exp,
X = df.cov,
iter = 10000,
verbose = FALSE,
varsel = TRUE,
family='gaussian')
# Extract pips
pips_ext_pfas <- ExtractPIPs(fit_ext_pfas)
pips_ext_pfas <- pips_ext_pfas %>% mutate(outcome = "Externalizing", group = "PFAS")
## PBDE
# Define dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi, `Substance Use`)
# Define mixture, covariates
df.exp <- mixture %>% dplyr::select(PBDE47.log, PBDE99.log, PBDE100.log)
df.cov <- mixture %>% dplyr::select(`Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi, `Substance Use`)
df.cov <- sapply(df.cov, as.numeric)
# BKMR
# Internalizing
set.seed(123)
fit_int_pbde <- kmbayes(y = mixture$Internalizing,
Z = df.exp,
X = df.cov,
iter = 10000,
verbose = FALSE,
varsel = TRUE,
family='gaussian')
# Extract pips
pips_int_pbde <- ExtractPIPs(fit_int_pbde)
pips_int_pbde <- pips_int_pbde %>% mutate(outcome = "Internalizing", group = "PBDE")
# Externalizing
set.seed(123)
fit_ext_pbde <- kmbayes(y = mixture$Externalizing,
Z = df.exp,
X = df.cov,
iter = 10000,
verbose = FALSE,
varsel = TRUE,
family='gaussian')
# Extract pips
pips_ext_pbde <- ExtractPIPs(fit_ext_pbde)
pips_ext_pbde <- pips_ext_pbde %>% mutate(outcome = "Externalizing", group = "PBDE")
Posterior inclusion probabilities (PIPs) reflecting the importance of each exposure in Bayesian Kernal Machine Regression models estimating the association between natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in first trimester serum samples, a mixture of PFAS alone, and a mixture of PBDEs alone.
Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.
# Combine results
bkmr <- rbind(pips_int, pips_ext, pips_int_pfas, pips_ext_pfas, pips_int_pbde, pips_ext_pbde)
# Table
table_s7 <- bkmr %>%
mutate_if(is.numeric, round, digits = 2) %>%
dplyr::select(outcome, group, variable, PIP) %>%
mutate(variable = recode(variable,
"PFHXS.log" = "PFHXS",
"PFOA.log" = "PFOA",
"PFNA.log" = "PFNA",
"PFOS.log" = "PFOS",
"PBDE47.log" = "BDE-47",
"PBDE99.log" = "BDE-99",
"PBDE100.log" = "BDE-100")) %>%
pivot_wider(names_from = outcome, values_from = PIP) %>%
as_grouped_data(groups = "group") %>%
flextable() %>%
add_body_row(values = c("", "", "PIP", "PIP")) %>%
set_header_labels(group = "", variable = "Chemical") %>%
font(fontname = "Tahoma") %>%
autofit() %>%
theme_alafoli()
table_s7
Chemical | Internalizing | Externalizing | |
|---|---|---|---|
PIP | PIP | ||
Total Sample | |||
PFHXS | 0.47 | 0.44 | |
PFOS | 0.27 | 0.31 | |
PFOA | 0.24 | 0.16 | |
PFNA | 0.34 | 0.32 | |
BDE-47 | 0.46 | 0.28 | |
BDE-99 | 0.37 | 0.27 | |
BDE-100 | 0.26 | 0.36 | |
PFAS | |||
PFHXS | 0.50 | 0.48 | |
PFOS | 0.22 | 0.32 | |
PFOA | 0.19 | 0.23 | |
PFNA | 0.35 | 0.41 | |
PBDE | |||
BDE-47 | 0.54 | 0.44 | |
BDE-99 | 0.56 | 0.46 | |
BDE-100 | 0.30 | 0.38 |
Cumulative effect (95% credible intervals) of the BKMR mixture of natural log transformed per- and polyfluoroalkyl substances (PFAS) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.
Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.
### TOTAL SAMPLE
df.exp <- mixture %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log)
# Internalizing
risk.overall_int <- OverallRiskSummaries(fit =fit_int, y = mixture$Internalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov),
qs = seq(0.25, 0.75, by = 0.05),
q.fixed = 0.5, method = "exact")
# Externalizing
risk.overall_ext <- OverallRiskSummaries(fit = fit_ext, y = mixture$Externalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov),
qs = seq(0.25, 0.75, by = 0.05),
q.fixed = 0.5, method = "exact")
# Merge
risk.overall_int <- risk.overall_int %>% mutate(outcome = "Internalizing")
risk.overall_ext <- risk.overall_ext %>% mutate(outcome = "Externalizing")
risk.overall <- rbind(risk.overall_int, risk.overall_ext)
## PFAS
df.exp <- mixture %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log)
# Internalizing
risk.overall_int_pfas <- OverallRiskSummaries(fit = fit_int_pfas, y = mixture$Internalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov),
qs = seq(0.25, 0.75, by = 0.05),
q.fixed = 0.5, method = "exact")
# Externalizing
risk.overall_ext_pfas <- OverallRiskSummaries(fit = fit_ext_pfas, y = mixture$Externalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov),
qs = seq(0.25, 0.75, by = 0.05),
q.fixed = 0.5, method = "exact")
# Merge
risk.overall_int_pfas <- risk.overall_int_pfas %>% mutate(outcome = "Internalizing")
risk.overall_ext_pfas <- risk.overall_ext_pfas %>% mutate(outcome = "Externalizing")
risk.overall_pfas <- rbind(risk.overall_int_pfas, risk.overall_ext_pfas)
## PBDE
df.exp <- mixture %>% dplyr::select(PBDE47.log, PBDE99.log, PBDE100.log)
# Internalizing
risk.overall_int_pbde <- OverallRiskSummaries(fit = fit_int_pbde, y = mixture$Internalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov),
qs = seq(0.25, 0.75, by = 0.05),
q.fixed = 0.5, method = "exact")
# Externalizing
risk.overall_ext_pbde <- OverallRiskSummaries(fit = fit_ext_pbde, y = mixture$Externalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov),
qs = seq(0.25, 0.75, by = 0.05),
q.fixed = 0.5, method = "exact")
# Merge
risk.overall_int_pbde <- risk.overall_int_pbde %>% mutate(outcome = "Internalizing")
risk.overall_ext_pbde <- risk.overall_ext_pbde %>% mutate(outcome = "Externalizing")
risk.overall_pbde <- rbind(risk.overall_int_pbde, risk.overall_ext_pbde)
# Plot
risk.overall_int_bkmr <- ggplot(risk.overall_int, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) +
geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey25") +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) +
ggtitle("Total Sample
Internalizing") +
theme_classic(base_family = "AppleGothic") +
labs(x = "Quantile", y = "Beta (95 % Credible Interval)") +
theme(text = element_text(size = 25, color = "grey25"),
axis.text = element_text(color = "grey25", size = 25),
legend.title = element_blank(),
plot.title = element_text(hjust = 0.5))
risk.overall_ext_bkmr <- ggplot(risk.overall_ext, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) +
geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey55") +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
scale_x_continuous(limits = c(0.25, 0.75), breaks = c(0.30, 0.40, 0.50, 0.60, 0.70)) +
scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) +
theme_classic() +
labs(x = "Quantile", y = "Beta (95 % Credible Interval)", title = "Externalizing") +
theme(text = element_text(size = 25, color = "grey25"),
axis.text = element_text(color = "grey25", size = 25),
legend.title = element_blank(),
plot.title = element_text(hjust = 0.5))
risk.overall_int_pfas_bkmr <- ggplot(risk.overall_int_pfas, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) +
geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey25") +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) +
ggtitle("PFAS
Internalizing") +
theme_classic(base_family = "AppleGothic") +
labs(x = "Quantile", y = "Beta (95 % Credible Interval)") +
theme(text = element_text(size = 25, color = "grey25"),
axis.text = element_text(color = "grey25", size = 25),
legend.title = element_blank(),
plot.title = element_text(hjust = 0.5))
risk.overall_ext_pfas_bkmr <- ggplot(risk.overall_ext_pfas, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) +
geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey55") +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) +
theme_classic(base_family = "AppleGothic") +
labs(x = "Quantile", y = "Beta (95 % Credible Interval)", title = "Externalizing") +
theme(text = element_text(size = 25, color = "grey25"),
axis.text = element_text(color = "grey25", size = 25),
legend.title = element_blank(),
plot.title = element_text(hjust = 0.5))
risk.overall_int_pbde_bkmr <- ggplot(risk.overall_int_pbde, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) +
geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey25") +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) +
ggtitle("PBDE
Internalizing") +
theme_classic(base_family = "AppleGothic") +
labs(x = "Quantile", y = "Beta (95 % Credible Interval)") +
theme(text = element_text(size = 25, color = "grey25"),
axis.text = element_text(color = "grey25", size = 25),
legend.title = element_blank(),
plot.title = element_text(hjust = 0.5))
risk.overall_ext_pbde_bkmr <- ggplot(risk.overall_ext_pbde, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) +
geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey55") +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) +
theme_classic(base_family = "AppleGothic") +
labs(x = "Quantile", y = "Beta (95 % Credible Interval)", title = "Externalizing") +
theme(text = element_text(size = 25, color = "grey25"),
axis.text = element_text(color = "grey25", size = 25),
legend.title = element_blank(),
plot.title = element_text(hjust = 0.5))
figure_4 <- plot_grid(risk.overall_int_bkmr, risk.overall_int_pfas_bkmr, risk.overall_int_pbde_bkmr, risk.overall_ext_bkmr,risk.overall_ext_pfas_bkmr, risk.overall_ext_pbde_bkmr, ncol = 3, row = 2)
figure_4
Univariate exposure-response functions indicating change in internalizing and externalizing behaviors resulting from individual natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in first trimester serum samples, while fixing remaining exposures in the mixture at the 50th percentile, estimated using Bayesian Kernal Machine Regression.
Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.
# Build labels
chemicals_label <- c("PFHXS", "PFOS", "PFOA", "PFNA", "BDE-47", "BDE-99", "BDE-100")
names(chemicals_label) <- c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")
# Internalizing
pred.resp.univar <- PredictorResponseUnivar(fit = fit_int)
int_univar <- ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) +
geom_smooth(stat = "identity", color ="grey25", size = 2) +
facet_wrap(~ variable, label = labeller(variable = chemicals_label)) +
scale_y_continuous(breaks = c(-1, 0, 1)) +
ylab("Internalizing") +
xlab("Chemical") +
theme_bw(base_family = "AppleGothic") +
theme(text = element_text(size = 30, color = "grey25"),
axis.text = element_text(color = "grey25", size = 30),
strip.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
# Externalizing
pred.resp.univar <- PredictorResponseUnivar(fit = fit_ext)
ext_univar <- ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) +
geom_smooth(stat = "identity", color = "grey25", size = 2) +
facet_wrap(~ variable, label = labeller(variable = chemicals_label)) +
scale_y_continuous(breaks = c(-1, 0, 1)) +
xlab("Chemical") +
ylab("Externalizing") +
theme_bw(base_family = "AppleGothic") +
theme(text = element_text(size = 30, color = "grey25"),
axis.text = element_text(color = "grey25", size = 30),
strip.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
figure_s2 <- plot_grid(int_univar, ext_univar, ncol = 2, align = "h")
figure_s2
Bivariate exposure-response functions indicating change in internalizing and externalizing behaviors resulting from individual natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in first trimester serum samples, while fixing remaining exposures in the mixture at specific percentiles, estimated using Bayesian Kernal Machine Regression.
Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.
# Define mixture
df.exp <- mixture %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log)
# Internalizing
pred.resp.bivar <- PredictorResponseBivar(fit = fit_int, min.plot.dist = 1)
pred.resp.bivar.levels <- PredictorResponseBivarLevels(
pred.resp.df = pred.resp.bivar,
Z = as.matrix(df.exp), qs = c(0.25, 0.5, 0.75))
# Plot
# Create labels
chemicals_label <- c("PFHXS", "PFOS", "PFOA", "PFNA", "BDE-47", "BDE-99", "BDE-100")
names(chemicals_label) <- c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")
pred.resp.bivar.levels$variable1 <- factor(pred.resp.bivar.levels$variable1, levels = c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))
pred.resp.bivar.levels$variable2 <- factor(pred.resp.bivar.levels$variable2, levels = c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))
figure_s21 <- ggplot(pred.resp.bivar.levels, aes(z1, est)) +
geom_smooth(aes(col = quantile), stat = "identity") +
facet_grid(variable2 ~ variable1, label = labeller(variable1 = chemicals_label, variable2 = chemicals_label)) +
scale_y_continuous(breaks = c(-0.4, 0, 0.4)) +
xlab("Chemical") +
ylab("Internalizing") +
scale_color_manual(values = c("steelblue3", "salmon3", "springgreen4")) +
theme_bw(base_family = "AppleGothic") +
theme(text = element_text(size = 30, color = "grey25"),
axis.text = element_text(color = "grey25", size = 30),
strip.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "bottom")
# Externalizing
pred.resp.bivar <- PredictorResponseBivar(fit = fit_ext, min.plot.dist = 1)
pred.resp.bivar.levels <- PredictorResponseBivarLevels(
pred.resp.df = pred.resp.bivar,
Z = as.matrix(df.exp), qs = c(0.25, 0.5, 0.75))
pred.resp.bivar.levels$variable1 <- factor(pred.resp.bivar.levels$variable1, levels = c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))
pred.resp.bivar.levels$variable2 <- factor(pred.resp.bivar.levels$variable2, levels = c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))
# Plot
figure_s22 <- ggplot(pred.resp.bivar.levels, aes(z1, est)) +
geom_smooth(aes(col = quantile), stat = "identity") +
facet_grid(variable2 ~ variable1, label = labeller(variable1 = chemicals_label, variable2 = chemicals_label)) +
scale_y_continuous(breaks = c(-0.4, 0, 0.4)) +
xlab("Chemical") +
ylab("Externalizing") +
scale_color_manual(values = c("steelblue3", "salmon3", "springgreen4")) +
theme_bw(base_family = "AppleGothic") +
theme(text = element_text(size = 30, color = "grey25"),
axis.text = element_text(color = "grey25", size = 30),
strip.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "bottom")
figure_s3 <- plot_grid(figure_s21,figure_s22, col = 2, align = "h")
figure_s3
Self-organizing maps (SOM) cluster star plot reflecting patterns of exposure to per- and polyfluoroalkyl substances (PFAS) and polybrominated diphenyl ethers (PBDEs), measured in serum samples obtained during the first trimester, and association between SOM clusters and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.
# Define mixture
df.exp <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log)
df.exp <- df.exp %>% rename(PFHXS = PFHXS.log, PFOS = PFOS.log, PFOA = PFOA.log, PFNA = PFNA.log, "BDE-47" = PBDE47.log, "BDE-99" = PBDE99.log, "BDE-100" = PBDE100.log)
# Explore clustering stats
# grp.str<-grpeval(df.exp, kmx=10)
# Initialize neurons/size of grid
grid <- somgrid(xdim = 3, ydim = 1, topo = "hexagonal")
# Train SOM model
set.seed(123)
som.mix <- som(X = as.matrix(df.exp), grid = grid, rlen = 500)
# Check if number of iterations is ok (should be stable)
# plot(som.mix, type = "changes")
## Visualize
# N
# plot(som.mix, type = "counts")
# plot(som.mix, type = "mapping", keepMargins = TRUE)
# Codes/weights
par(cex = 8)
plot(som.mix, type = "codes", main = "SOM Clusters")
# Get unit classification
mix$unit_clusters <- som.mix$unit.classif
# get events (weight vector)
# getCodes(som.mix)
# Regression based on unit classifications
# Assign reference group
mix$unit_clusters <- relevel(as.factor(mix$unit_clusters), ref = "1")
# Internalizing
som_int_reg <- lm(Internalizing ~ unit_clusters + `Mean Child Age` + `Infant Sex` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`, data = mix)
som_int_reg_df <- tidy(som_int_reg , conf.int = TRUE) %>% filter(term %in% c("unit_clusters2", "unit_clusters3")) %>% dplyr::select(term, estimate, conf.low, conf.high) %>% mutate(outcome = "Internalizing")
# Externalizing
som_ext_reg <- lm(Externalizing ~ unit_clusters + `Mean Child Age` + `Infant Sex` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`, data = mix)
som_ext_reg_df <- tidy(som_ext_reg , conf.int = TRUE) %>% filter(term %in% c("unit_clusters2", "unit_clusters3")) %>% dplyr::select(term, estimate, conf.low, conf.high) %>% mutate (outcome = "Externalizing")
som_reg_df <- rbind(som_int_reg_df, som_ext_reg_df )
som_reg_df$outcome <- factor(som_reg_df$outcome, levels = c("Internalizing", "Externalizing"))
figure_5 <- ggplot(som_reg_df, aes(x = estimate, y = term, color = outcome)) +
geom_pointrange(aes(xmin = conf.low, xmax = conf.high), size = 4, linewidth = 4) +
geom_vline(xintercept = 0, color = "grey25", linetype = "dashed", alpha = 0.5) +
facet_grid(~outcome, scales = "free") +
xlab("Beta (95% Confidence Interval)") +
ylab("Chemical") +
scale_x_continuous(n.breaks=4) +
scale_y_discrete(labels = c("Cluster 2", "Cluster 3")) +
scale_color_manual(values = c("grey25", "grey55")) +
coord_flip() +
theme_classic(base_family = "AppleGothic") +
theme(title = element_text(size = 50, color = "grey25"),
axis.title.x = element_text(size = 50, color = "grey25"),
axis.title.y = element_text(size = 50, color = "grey25"),
axis.text.x = element_text(size = 50, color = "grey25", vjust = 0.6),
axis.text.y = element_text(size = 50, color = "grey25"),
strip.background = element_rect(fill="white"),
strip.text = element_text(colour = 'grey25', face = "bold", size = 50),
legend.position = "none",
strip.background.x = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank())
figure_5
Linear regression associations between self-organizing map (SOM) clusters and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.
Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.
# Table
table_s8 <- som_reg_df %>%
mutate_if(is.numeric, round, digits = 2) %>%
mutate(Beta_CI = paste0(estimate, " (", conf.low, ", ", conf.high, ")")) %>%
mutate(term = recode(term,
"unit_clusters2" = "Cluster 2",
"unit_clusters3" = "Cluster 3")) %>%
dplyr::select(term, outcome, Beta_CI) %>%
pivot_wider(names_from = outcome, values_from = (Beta_CI)) %>%
flextable() %>%
add_body_row(values = c("" ,"Beta (95% CI)", "Beta (95% CI)")) %>%
font(fontname = "Tahoma") %>%
set_header_labels(term = "", sex = "Sex") %>%
autofit() %>%
theme_alafoli()
table_s8
Internalizing | Externalizing | |
|---|---|---|
Beta (95% CI) | Beta (95% CI) | |
Cluster 2 | 0.41 (-0.01, 0.83) | 0.5 (0.09, 0.91) |
Cluster 3 | 0.44 (0.08, 0.81) | 0.54 (0.18, 0.89) |